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A linear-time algorithm is presented for the construction of the Gibbs distribution of configu- 
rations in the Ising model, on a quantum computer. The algorithm is designed so that each run 
provides one configuration with a quantum probability equal to the corresponding thermodynamic 
weight. The partition function is thus approximated efficiently. The algorithm neither suffers from 
critical slowing down, nor gets stuck in local minima. The algorithm can be applied in any dimen- 
sion, to a class of spin-glass Ising models with a finite portion of frustrated plaquettes, diluted Ising 
models, and models with a magnetic field. 



I. INTRODUCTION 

The algorithm of Shor jlj] for the polynomial time so- 
lution of the factorization problem on a quantum com- 
puter was received with much excitement in the com- 
puter science and physics communities ^,|) . It indicates 
that quantum computers have a potential for the effective 
solution of problems which are unfeasible on a classical 
computer. The actual utilization of this potential would 
require, in addition to many years of work on the hard- 
ware, the development of algorithms which would opti- 
mally exploit the strengths while overcoming the short- 
comings of quantum computers, in particular the prob- 
lem of decoherence ^|||- Unlike classical computers in 
which each bit is a two state system which can be either 
in state or 1 the quantum bit (or qubit) can be in any 
superposition of the form: 



|^) = a |0)+ai|l) 



(1.1) 



as long as | cko | 2 + |c*i| 2 = 1- When a measurement of 
the qubit takes place the result will be the state |0) with 
probability |ao| 2 or |1) with probability |ai| 2 . In the first 
case the system will then remain in the state |0) while 
in the second case it will remain in the state Due 
to superposition, a system of N qubits is described by 
a unit vector in a 2^ dimensional complex vector space 
(the Hilbert space) of the form: 



N 



E 



(1.2) 



where \i) are the 2^ basis vectors and \ a i\ 2 = 1- Com- 
putations are done by changing the state of the system. 
Since conservation of probability is required only unitary 
transformations are allowed. One source of the poten- 
tial strength of quantum computers is due to the fact 



that computations are performed simultaneously on all 
2 N states in the superposition. This amounts to an ex- 
ponential parallelism compared to classical computers in 
which at any given time only one number can be stored 
in the ./V bit register. However, this feature alone is not 
enough since only a single quantum state is available after 
measurement ||. The full power of quantum computa- 
tion is realized only when the superposition principle is 
combined with the ability of quantum states to interfere. 
The latter has no classical analogue, and is the quantum 
ingredient which allows one to selectively control which 
state will have the highest probability of appearing after 
measurement. 

Similarly to classical computers, it turns out that all 
unitary transformations involving N qubits can be bro- 
ken into two qubit unitary transformations j?] ^| . This 
allows one to construct a universal set of binary gates 
which is capable of implementing all the required oper- 
ations. The actual construction of a quantum computer 
is a formidable task that is believed to require many 
years before a basic prototype will be ready. Some of 
the potential physical media proposed for quantum com- 
puters include ions in ion traps jOj, atoms coupled to 
optical resonators [fl2|| , interacting electrons in quantum 
dots |L3]], and Ramsey atomic interferometry |l4[ . 

The main difficulty identified so far in the construc- 
tion of a quantum computer is the decoherence of the 
quantum superposition due to interaction with the envi- 
ronment. To avoid errors one needs to isolate the quan- 
tum computer from the environment as much as possible. 
Some redundancy combined with error correcting codes 
is considered as a promising way to reduce the accumu- 
lation of errors during the computation jl5]-^2| . Another 
potential difficulty is to maintain sufficient precision so 
that the quantum computer will provide accurate results 
even after many steps of computing. Therefore, an ef- 
ficient quantum algorithm should satisfy not only that 
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the time and memory required for the computation grow 
polynomially with the input size, but also the precision: 
the number of bits of precision should grow only loga- 
rithmically in the input size Jlp^| . 

In addition to the prospects of building a quantum 
computer, the experimental work stimulated by this field 
is expected to provide new insights on the foundations 
of quantum mechanics, as well as to lead to progress 
in the development of new technologies. Furthermore, 
the new perspective that quantum computers provide 
about the complexity of algorithms is highly valuable 
for the theory of computation even without the physical 
construction of such a computer. In particular, Shor's 
algorithm for the polynomial time solution of the fac- 
torization problem Q| has shown that problems which 
are considered intractable on classical computers may 
be tractable on a quantum computer, although some re- 
strictions are known related to the potential power of 
a quantum computer p3[ |. Feynman pi] was the first 
to suggest that quantum computers might be exponen- 
tially faster than classical computers at simulating quan- 
tum mechanical systems with short range interactions. A 
general demonstration to this effect was given by Lloyd 
P5[ , who also argued that quantum computers could ef- 
ficiently calculate spin-spin correlation functions in Ising 
models. Some explicit algorithms were later proposed 
for simulating physical systems on quantum computers. 
These include the Schrodinger equation for interacting 
many-body systems [p6|- p9|, t he Dirac equation |}0| and 
the quantum baker's map |3l| . 

In this paper we consider a broad class of statistical 
physics problems involving Ising spin systems. We de- 
velop an algorithm for simulating these systems on a 
quantum computer. Here we define the Ising spin sys- 
tems and briefly review the numerical techniques in use 
for their simulation on classical computers. These sys- 
tems are important both as models of magnetic phase 
transitions and as the most useful models for the analyt- 
ical and numerical studies of phase transitions in general 
p2p^ ]. The numerical simulation of such systems has 
been an active field of research for the last five decades 
since the introduction of the Metropolis algorithm |34[ ]. 
Typically, spin systems are defined on a d dimensional 
lattice in which there are N spins s^, one at each lat- 
tice site i, and nearest neighbor coupling between spins. 
The energy of the system is given by the nearest-neighbor 
Edwards- Anderson Hamiltonian (3^] : 

N 

H = - } JjjSjSj - } hjSj, (1.3) 
(i,i> *=1 

where (i, j) represents summation only over pairs of near- 
est neighbors, Jy is the coupling between i and j, and 
hi is the local magnetic field. The most commonly stud- 
ied model is the Ising model |36| in which each spin has 
two states Si = ±1. The bonds then determine the 
nature of the interactions. In the ordinary ferromagnetic 



(antiferromagnetic) Ising model all bonds satisfy J L j = J 
(Jij = —J) where J > 0. In the ±J Ising spin glass there 
are quenched random bonds chosen from a bimodal dis- 
tribution PiJij) =p5(Jij — J) + (1 — p)S(J ij + J). The 
random bonds in the Ising spin glass can also be drawn 
from a continuous distribution such as the Gaussian dis- 
tribution. 

Numerical studies of spin systems have been performed 
for a vast variety of lattices including the square and tri- 
angular lattices in two dimensions (2D), cubic, hexag- 
onal and hexagonal closed packed in three dimensions 
(3D). Here we will concentrate mainly on finite hypercu- 
bic lattices in d dimensions, which include N — L d sites, 
where L is the number of sites along each edge. Since 
each spin has two states, these Ising spin systems exhibit 
an exponentially large phase space of 2 N configurations. 
The partition function of Ising type models is given by 
a sum over all configurations: Z = Y^{ s } ex P(~/^[{ s }]) 
where (3 = 1/(A:bT), ks is the Boltzmann constant and 
T is the temperature. For a system in thermodynamic 
equilibrium, the partition function provides the statis- 
tical weight of each one of the 2 N configurations. The 
statistical weight of the configuration {si},i = l, . . . , N 
is given by: Pr[{si}] = exp(— (3H[{si}])/Z. Therefore, 
if the partition function is known one can obtain exact 
results for all the thermodynamic quantities such as the 
magnetization, susceptibility and specific heat. Models 
for which analytical calculations of this type can be per- 
formed include a variety of one dimensional (ID) models 
and the 2D Ising model However, for most systems 
of interest, including the 3-D Ising model and most Ising 
spin glass models, no analytical calculation of the parti- 
tion function is available Pq] . 

The size of the input in computations of Ising spin sys- 
tems is simply the number of spins N plus the number of 
bonds Nb which connect these spins. In models of short 
range interaction considered here the number of bonds 
is of order O(N). An exact numerical calculation of the 
partition function or any thermodynamic quantity would 
involve a summation over the 2 N terms which appear 
in the partition function. As the system size increases 
the computation time would increase exponentially, and 
this is obviously not feasible. In order to obtain ther- 
modynamic averages a variety of Monte Carlo methods 
have been developed. These methods involve a sequen- 
tial random sampling of the phase space moving from 
one configuration to the next according to a properly de- 
signed Markov process. In order to sample all configura- 
tions with the appropriate thermodynamic weights, the 
Markov process must be able to access the entire phase 
space and to satisfy the detailed balance condition: 

Pr({ Sl }) • W({ Sl } - {*}') = Pr({ Sl }') • W({ Sl }' - { Si }) 

(1.4) 

where {sj} and {si}' are any two states of the system 
and PF({si} — » {si}') is the transition probability from 
one state to the other in a single move of the Markov 
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process |3^] . When these conditions are satisfied one can 
use the Monte Carlo results to obtain approximations to 
thermodynamic quantities. 

In the most commonly used Metropolis algorithm [pJJ , 
at each time step one spin is chosen randomly. The en- 
ergy difference AE which would occur due to flipping 
the chosen spin is calculated. If AE < the move is 
accepted and the spin is flipped. If AE > the move 
is accepted with probability p — e~@ AE . Since this rule 
satisfies detailed balance, one can take samples of the 
configurations during the run to obtain properties of the 
equilibrium phases such as magnetization, susceptibility, 
correlation function, correlation length, etc. Although 
for large systems it is feasible to scan only an exponen- 
tially small part of the phase space, this part typically 
has an exponentially large weight and therefore, Monte 
Carlo simulations provide good approximations for the 
thermodynamic quantities. 

The Metropolis algorithm and related techniques 
which involve flipping one spin at a time are efficient 
as long as the system is not too close to a critical point, 
i.e., a second order phase transition. Near the critical 
point the simulation suffers from "critical slowing down" 
and the number of Monte Carlo steps needed between 
uncorrelated configurations grows as L z where z is the 
dynamical critical exponent |B3| . The reason for this 
slow dynamics is that near the critical point there are 
large clusters of highly correlated spins. In this situa- 
tion there is a very small likelihood of flipping an entire 
cluster due to the large energy barrier involved. To over- 
come this difficulty, cluster algorithms were introduced, 
in which entire clusters are flipped at once, in a way that 
maintains detailed balance p9[ |. 

In addition to the regular lattice spin systems, there 
has been much interest in the study of disordered systems 
such as frustrated Ising models JlO 41 1 and the Ising spin 
glass ^,ff3| • These systems exhibit competing ferromag- 
netic and antiferromagnetic interactions. In particular, 
in plaquettes which include an odd number of antiferro- 
magnetic bonds it is not possible to satisfy all the bonds 
simultaneously and the system is thus frustrated |io| , pl| . 
Spin glasses exhibit a complex energy landscape with a 
large number of metastable states or local minima. Since 
these minima are separated by energy barriers, when the 
system is simulated using Monte Carlo methods at low 
temperatures, it tends to be trapped around one local 
minimum from which it cannot escape. When the simu- 
lation is done at high temperature, the system can easily 
switch from the vicinity of one local minimum to an- 
other but cannot resolve the details, namely reach the 
local minimum itself. The approach of simulated anneal- 
ing pil in which the temperature is repeatedly raised 
and then slowly decreased was found useful for obtain- 
ing thermodynamic averages. In particular, it provides 
a probabilistic algorithm for exploring the local minima 
and for searching for the ground state of the system. The 
problem of finding the ground state of the short range 3D 
Ising spin glass, as well as the fully antiferromagnetic 2D 



Ising model in the presence of a constant magnetic field 
was shown by Barahona to belong to the class of non- 
deterministic polynomial time (NP)-hard problems, by a 
mapping to problems in graph theory |^5|-[48|. 

In this paper we present an algorithm for the study 
of a class of random-bond Ising spin systems on a quan- 
tum computer. By use of interference, the algorithm can 
construct, with linear complexity, a lattice with a fixed 
portion of plaquettes having predetermined bonds. The 
bonds connecting plaquettes are determined randomly, 
with the probability of obtaining a non-frustrated inter- 
mediate plaquette being higher than that of a frustrated 
one. The superposition property of a quantum com- 
puter can be used in order to include the entire phase 
space of the resulting Ising system, such that the quan- 
tum mechanical probability of each one of the 2 N states 
equals the thermodynamic weight of the corresponding 
spin configuration. In this sense the algorithm is exact. 
Once such a superposition is constructed one can perform 
a measurement of all spins, which provides one of the 
2 N configurations. Since the probabilities of the quan- 
tum states are ordered by the thermodynamic weights of 
the corresponding spin configurations, the partition func- 
tion is constructed efficiently. Putting aside questions of 
degeneracy, the lower the energy of the configuration, 
the more likely it is to be obtained upon measurement. 
Therefore the algorithm can be used to find ground state 
configurations of the spin system. Unlike Monte Carlo 
simulations which require a minimal number of steps be- 
tween measurements to reduce the autocorrelation func- 
tion between consecutively measured configurations |53| , 
measurements on the quantum computer are totally un- 
correlated since the superposition is constructed anew 
before every measurement. Successive runs of the algo- 
rithm therefore involve no dynamics of moving from one 
configuration to another. A situation in which this pro- 
cedure is especially useful, is in the vicinity of a critical 
point, where Monte Carlo simulations may suffer from 
critical slowing down. While cluster algorithms have 
been invented for regular Ising models [p9|j4g| l which es- 
sentially solve this problem, they are very limited in scope 
and can treat essentially only Ising systems with periodic 
bond-structure. The present algorithm is more general in 
that it applies to random-bond Ising models as well, and 
avoids the problem of critical slowing down altogether. 

The paper is organized as follows. The construction 
of the superposition of states for the ID Ising model, 
with quantum probabilities equal to the thermodynamic 
weights is considered in Section ||. Higher dimensional 
Ising systems including the Bethe lattice, the 2D and 3D 
Ising models are conside red in Section 



III. A magnetic 



field is introduced in Sec, IV, The conclusions appear in 
Section [v|. 
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II. ID ISING MODEL 

We begin our exposition of the algorithm by treat- 
ing the simple case of a ID Ising model. Starting from 
the fully ferromagnetic open chain, we will gradually in- 
troduce complexity, by considering the antiferromagnetic 
case, mixed ferro-antifcrromagnctic models, spin-glasses, 
and finally close the boundary conditions. This last op- 
eration, which enables the use of transfer matrices in the 
classical case, allows for a comprehensive treatment of 
the 2D and higher dimensional models. 



A. The Ferromagnetic Case 

The Hamiltonian for a linear, open ferromagnetic sys- 
tem of N spins Si = ±1 is [ jspf : 



N-l 



(2.1) 



i=i 



where J > 0. Let y e [0, 2 N -l], and {s}% be the iV-digit 
binary expansion of y using Sj = —1 for and Si = +1 
for 1. The notation {s}y can also denote one of the 2 
spin configurations, with thermodynamic weight: 



?r[{s} y ] = 4+ 



-f3H+[{s} y ] 



Z 



(2.2) 



N 



(the N superscript on {s}y will be suppressed where it 
is obvious) where: 



7+ 

N 



E 



(2.3) 



is the partition function. It is the task of the algorithm 
to exactly calculate the probabilities above, in a manner 
which allows an easy identification of the configuration 
whose probability was found. To this end we introduce 
an iV-qubit register {s}y — \si, S2, sjv), where now 
Si = ±1 denote the ground and excited states of the i th 
qubit (it is convenient to use the same notation for the 
classical spins and the qubits, and this should not cause 
any confusion). We term this a register of "spin-qubits" . 
Let |{— }n) denote the ground state of all spin-qubits. 
We seek a unitary operator which simulates the Ising 
model in the following sense: 



\(My\n\i 



}n)\ =Pr[{*U 



(2.4) 



Thus Tjy evolves the qubit register into a superposition 
in which every state uniquely codes for an Ising configu- 
ration of spins, with a quantum probability equal to the 
thermodynamic weight of that configuration. must 
be a "valid" quantum computer operator, i.e., it must 
be decomposable into a product of a polynomial (in N) 
number of 1- and 2-qubit unitary operators only 0] . Such 



a decomposition is possible with the following two oper- 
ators: a 1-qubit 7r/2 rotation, 



Ri\si, 
1 „ 
V2 

and a 2-qubit "Ising-entanglement" : 

Sij\ s li Si, Sj, Sjv) = 

1 

. , Si , . . . , Sj , . . . , Sjv 



., Si, Sn) — 

|si, ...,-Si, ...,sjv) 



Si\si, Si, sn)), (2.5) 



-M^lsx, 



JSi 



S\ , . . . , Si 



, sjv)) , 



where 



cee-Z+ = 2cosh(/3J). 



(2.6) 



(2.7) 



In what follows we will suppress the full register and 
indicate only the qubits operated on. It is straightfor- 
ward to check that Ri and Sfj are indeed unitary, e.g., 
by considering their matrix representations in the ba- 



sis where |-) = (1,0), |+) = (0,1), | - -) 
|-+) = (0,1,0,0),|+-) = (0,0,1,0),|++) 



(1,0,0,0), 
(0,0,0,1): 



R 



and 



V 











V2 ( 









x J 


x J 








—x~ J 


X J 














x~ J 


x J 





- 


-x J 


X~ J 



(2.9) 



It is interesting to note the similarity to the classical ID 
transfer matrix, 



,.2J 



„-2J 



(2.10) 



The x J vs x 2J comes from the fact that in S + we have 
amplitudes, not probabilities. 

The operator which simulates the ID Ising problem 
can now be written as: 



1 N 



=iV-l 



(2.11) 



Thus Tjy is a ir/2 rotation of the first qubit, followed by 
Ising entanglements of successive pairs of qubits. This 
bares some resemblance to the procedure using a trans- 
fer matrix to solve the ID Ising model. The num- 
ber of required operations is exactly N. The general 
"recipe" for writing down this operator (in the absence 
of closed loops) is the following: one always applies a 
7r/2 rotation to the first qubit, and then substitutes an 
Ising-entanglement operator for each interacting pair of 
nearest-neighbor spins in the Hamiltonian. 
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3 


4 


5 


6 


Sl,S 2 , S 3 , S 4 ) 


initial 


Ri 


°12 


°23 


S + 

°34 


\ 


1 


1/V2 


X J /V2c 




x 3J /V2c 3/2 


h> 














-x'7V2c 3/2 


1" -> 











-1/V2c 


-x-i/V2c 3 ' 2 


h +) 














x J /V2c^ 2 


1 } 








-a;" J /^2c 


-a;- 27 /V2c 


-x- J /V2c i/2 


— H h) 














x~' iJ /^2c A ' 2 


- + + — > 











1/V2c 


x- J /V2c^ 2 


- + + +} 














-x J /V2c a/2 


H ) 





1/V2 


x~ J /V2c 


1/V2c 


x J /V2c^ 2 


H h) 














-x-'l^' 2 


1 H H-) 











-aT 27 /^ 


-x- 3J /V2c^ 2 


1 H H +> 














x- J /V2c 3 ' 2 


+ H ) 








-x J /V2c 


-l/\/2c 


-x J /V2c^ 2 


+ H h) 














x- J /V2c 3 ' 2 


+ + + -) 











x 27 /^ 


x J /V2c^ 2 


+ + + +) 














-x :iJ /V2c :i ' 2 



TABLE I. Amplitudes of register states for ID Ising model with up to 4 spins. 



It might be helpful to give an example at this point. For 
an open chain of N < 4 spins, Table || gives the ampli- 
tudes of four spins, at e ach stage of the algorithm, as 
calculated from Eq.( 2.11 ). It is easily verified that the 
squares of the amplitudes given in columns 4, 5, 6, agree 
with the thermodynamic weights [given by Eq.(2.1)] for 
N=2, 3, 4 spins respectively, with a ferromagnetic in- 
teraction. (In order to check, e.g, for N=2, ignore the 
entries for S3 and S4.) 



N 



We proceed to prove that Tjy indeed satisfies Eq.(2.4), 



by mathematical induction. Let us consider first the min- 
imal case N — 2 (for N = 1 we cannot apply S + ). We 
have T2 = S^-Ri, so that 

n\{-h) = ^su\ -->+!+-» = 



— r 



-s 2 x JsiS2 \sus 2 ) 



-}-x J \ + +))} = 

(2.12) 



On the other hand, according to Eq.(2.2) the thermo- 
dynamic weights of these four states are, respectively, 
ef J /Z+, e-* 3 - 1 /Z+, e-* 3 - 1 /Z+, e f3J /Z+. These are ex- 
actly the squares of the above amplitudes. Next, assume 
by induction that Eq.(|]J) holds for Tj^, i.e., that 



where: 



2 N -1 

n\{-} N ) = e 4K s w. 



1 -H+{{S} V ], 



(2.13) 



>N 



(2.14) 



Using this, we must show that T^ +1 satisfies Eq.(2.4) for 
the ID Ising model with N + 1 spins. Now, 

t n+i\{-}n+i) = S NjN+1 T N \{-} N )\-) = 
2"-l 

E A l S N,N+l\{s}y)\-)> (2-15) 

y=o 

where the last equality follows from the induction hy- 
pothesis. But by definition of S + , 

$U + lK»}?>|-> = \ (x- ISN \{s}y)H ~ ^M^)) 

Inserting this into Eq.( 2.15| ), we have: 
7jV+il{-}jv+i) = 



E 



r -0H+Us\„]/2 1 



»=0 JZ 



s N f}J/2\ UX N 



|{<,+>). (2.16) 



The two terms in this sum arise from the two pos- 
sible states of |sjv+i), so writ ing the exponents as 
exp(sNSN+iPJ/'2) and using Eq.(2.1), we obtain: 



T+ +1 |{-} N+1 ) = £ -^e-f^^\{ S }^). 

y=0 JZ N + 1 



(i 
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The state \{s}£ +1 ) thus appears with a probability equal 
to its thermodynamic weight, which completes our proof. 
A useful corollary is that 

Corollary 1 Upon a measurement following the execu- 
tion of the algorithm, a state appears with a probability 
equal to the thermodynamic weight of the corresponding 
spin configuration. 

This implies that the present algorithm provides an 
approximation to the partition function which converges 
rapidly in the number of measurements. 



T n+1 \{-}n+i) - 

2 N -1 I— 

/ r / — V c 

y=o y zj n 

(* SN \U}y,-)-X- SN \My,+)) = 
2 jv + i_ 1 

£ ^L e -^ +1 [W,]/2| {s} ^+i). (2.21) 

y=0 ^JZ N + 1 

This proves that simulates the ID antiferromagnetic 
Ising model. 



B. The Antiferromagnetic Case 

In the antiferromagnetic case, the Hamiltonian for a 
linear, open system of N spins is 



JV-l 



H N — J s^Si+i, 



(2.17) 



where J > 0. To treat this case we define a properly 
modified version of the Ising-entanglement operator S + : 



S t j \si,sj 



1 



) = —j= (x""'\S i ,S : j) -r s;.r 
v C 



-J 8, 



K-Sj)), (2.18) 



In matrix form, 



/ X 3 X J 
-x J X~ J 




\ 


x J x~ J 



(2.19) 



\0 -x- J x J J 
and it is easily checked that S~j is unitary. We claim that 



-N 



-N-l 



Ri 



(2.20) 



simulates the antiferromagnetic "ID Ising model" , in the 
sense of Eq. (^J) . The proof is entirely analogous to the 
ferromagnetic case. First one considers the case N = 2, 
for which: 



1 



T 2 -\{-h) = ^=Su(| 



-)) 



[(x- J \ --)- x J \ +)) + (i^ + -) - x- J \ + +))] . 



The state-probabilities are just the ther modynamic 
weights which can be obtained from Eq.( [2.17 ) for N = 2. 
Repeating the induction argument that led to Eq.(2.16) 
one finds here: 



C. The "ID Spin-Glass" Case 

We now consider the simplest case of a random- 
bond Ising model with open boundary conditions: the 
quenched, mixed ferro-antiferromagnetic linear chain 
(also known as the ±J spin glass). The Hamiltonian 
in this case may be written as: 



N-l 



H 



N 



= -£•*■ 



(2.22) 



where J = (Ji, J2, J N-l) is a fixed set of parameters, 
each of which can be ±J and thus determines whether 
the interaction between Si and Sj+i is ferro- (+) or anti- 
ferromagnetic (— ). There are a total of 2 N ~ 1 J's for the 
length- N Ising chain, each of which can be regarded as 
a different realization of quenched disorder. The opera- 
tor which simulates the corresponding Ising problem is a 
natural generalization of T^: 



T 3 N 



n 5 

i=N-i 



Ji 



Ri, 



(2.23) 



where 151 



-JiSi 



I Sj x I S 1 ^ ' 



(2.24) 

The unitarity of S^f follows from the unitarity of 
and S~j. Let us prove that T N simulates the appropri- 
ate Ising problem, again, by induction: for N = 2 there 
are two realizations of the quenched disorder: J\ = ±J. 
Accordingly, there are two T-f's: T 2 + and T 2 ~. We have al- 
ready shown that these operators solve their correspond- 
ing Ising problem. Assume by induction that T N simu- 
lates the Ising problem for N spins. There are now four 
possibilities in going to N + 1, since both Jn-i and Jn 
can be ±J. In fact we have already dealt with the two 
cases Jn-i = Jn hi proving the algorithm for the fully 
ferro- and antiferromagnetic cases. But instead of con- 
sidering separately the cases Jn-i 7^ Jn, it will be more 
convenient to proceed generally. From the induction hy- 
pothesis we have: 
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2 N -1 

t j n \{-}n) = £ 

v=o 

A J = <t>N x -Hl\U}J_ 



(2.25) 



= Sl N AT ,Ti so that using Eq.fl2.24Q 



Now, T N+1 — ' \ _ \ | | - x - 



T 3 n+1 \{-}n)H 



2"-l 

E 

y=o 



N x -X J A{s}y] _L x 

J \/c 



(*- J ™lto^->-z Jw " w IW?,+>) = 



2 « -!-*_! 

E 



0AT+1 p -/?g^ + 1 [{ S } H ]/2|{ s }jVH 



7 J 



(2.26) 



The amplitude-squared of the configuration coded by 
{s}y +1 is exactly its thermodynamic weight for a given 
quenched disorder J, so this completes the proof for the 
mixed ferro-antiferromagnetic case. Of course, the fully 
ferro- and antiferromagnetic cases are specific instances 
of this general model. The implementation of t he algo- 
rithm in the present case, according to Eq.( [2.23 ), would 
entail using (apart from the ir/2 rotation), two differ- 
ent operators, in an order dictated by the sequence of 
ferro/antiferromagnetic bonds in the Ising model one 
wishes to solve. The complexity, however, remains O(N). 

The generalization to continuous interactions is 
straightforward. In this case: 



N-l 

Hn = — y GiSiS 



iSiSi+x , G = (Gx, G2, Gjf-i), 

(2.27) 

where now G is a set of independent random variables, 
which do not have to be identically distributed. Suppose 
one prepares a set G. This corresponds to choosing a 
certain realization of quenched disorder in the Ising spin 
glass. The quantum operator which simulates for the 
thermodynamic weights in this case is 



n * 

.i=N-l 



Gi 



(2.28) 





(x G < 










1 


-x~ G ' 


x Gi 








y/Ci 








x~ G * 


X 









-x Gi 


X 



-Gi 



(2.30) 



The only difference from of the ± J spin glass is that 
each Ising-entanglement now has its own normalization 
factor, which clearly has no effect on the proof of the 
algorithm. 

As iV increases finite size effects diminish. By con- 
struction, our algorithm will provide in 0(N) steps a 
ground state {s}* of a "ID spin glass" with probability: 



p* = exp(-/3H%l{s}*})/Z% 



7G 



(2.31) 



where: 



This probability can be made arbitrarily close to 1 by 
performing the construction for low enough tempera- 
ture. If a local minimum is found instead of the ground 
state, the entire process should be repeated until the 
ground state is obtained. What is the average number 
of steps required for locating the global minimum in this 
manner? After the first run one has probability p* of 
having found {s}*. If not, one failed with probability 
q = 1 — p* and then succeeded with probability p*, etc. 
Clearly the resulting distribution is geometric, and thus 
the average number of runs until the global minimum is 
found, is (n) = 1/p* 152). The total number of steps 
is seen to be 0(N)/p*. The question of the complexity 
of the algorithm for locating a ground state thus boils 
down to the scaling of p* with N . In ID this is in fact 
trivial: there are exactly two degenerate ground states, 
related by s% — > — s, Vi, obtained by simply following 
along the chain and satisfying all bonds. Their energy is 
E = —Nb J (where Nb = N — 1 is the number of bonds) 
so p* = exp(/3iVs J)/Zjv. The temperature appears here 
as a control parameter: let A be the difference in poten- 
tial energy between a ground state and the next lowest 
states {s}^. Then Eq.( [2.31| ) predicts that p* will become 
dominant, since p* /p' = cxp[A/(fcsT)]. This indicates 
that the probability of obtaining a ground state can be 
made arbitrarily close to 1 |53|. However, this is only 
true as long as the degeneracy g ms of metastable states 
with energy close (of the order of the average interaction 
strength (G)) to the ground state energy remains small in 
some proper sense. For higher dimensional spin glasses, 
it is well known that this number is iV-dependent. Addi- 
tional complications may arise in relation to the connec- 
tivity of the system in higher dimensions. We will return 
to this issue in later sections. In order to be able the dis- 
cuss such systems, the problem of closing the boundary 
conditions must first be discussed. This is done next. 



•Sr.' |.' 



Gi 
ij 
1 

Ci = 2cosh(/3Gi). 



(x- G *'*\8 i ,8 j ) + 8jX G » t \8 i ,-8 j )) 



In matrix form: 



D. Closing the Boundary Conditions 

It should be remarked that the algorithm as described 
(2.29) so far is in fact classical. A classical probabilistic com- 
puter can run the algorithm with exactly the same effi- 
ciency simply by randomly choosing spin to be up or 
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down with a probability determined by the thermody- 
namic weight of the configuration of all other i — 1 spins 
and bonds. The difference is of course that the clas- 
sical computer cannot store all 2 N spin configurations. 
However, this by itself does not enhance the computing 
power since only one configuration is accessible by mea- 
surement of the quantum register. The effective classical- 
ity of the algorithm is due to the fact that so far we have 
only employed superpositions. In Sec. II F we will employ 
the purely quantum effect of interference in order to deal 
with the problem of a ID Ising chain with closed bound- 
ary conditions. Here we introduce the operator required 
for closing a ID chain. Such a chain has as Hamiltonian 
in the ± J spin glass case: 



N 



i=l 

sn+i = si, J = (Ji, J2, ■ Jn)- 



(2.32) 



A reasonable approach to closing the loop on a quantum 
computer might seem to be an application of the opera- 
tor S N N 1 after T N . However, this does not work since it 
changes the amplitude of \si) which was already the cor- 
rect thermodynamic weight. It turns out that a different 
approach is needed. Instead of working on \si), one has to 
first introduce a work-qubit, say \w), on which the Ising- 
entanglement operation is performed: S^ 1 ^ ' . This places 
| to) in a superposition of up and down states. Closing the 
loop is then performed by comparing the state of \s\) to 
that of the workbit (rather than to that of |sjv)), which 
acts effectively as the fictitious spin sjv+i- If to = Si, 
the loop is closed, since this corresponds to sn + ± = si, 
as should be. However, one is then left to wonder what 
to do if w = —8%. Instead of discarding this possibility 
as improper, it turns out to be fruitful to adopt a more 



general point of view. As will be shown in Eq.(2.38), the 
following interpretation also holds: If to = si, the loop 
is closed ferromagnetically [sign(Jjy) = 1]; if to = — Si, 
the loop is closed antiferromagnetically [sign(JAr) = — 1]. 



Since the sign of the interaction is determined randomly, 
we can only specify the absolute value of Jn, hence the 



notation S 



\Jn\ 

N,w ' 



The comparison operation can be per- 



formed by an exclusive-or (XOR): 

X ij | Si , Sj ) — SiSj | Si , SiSj ) , 

or, in matrix form: 



X = 



0-1 
10 
0-10 
1 



(2.33) 



(2.34) 



Thus following S^Z 1 one applies Xx, w \ the combined op- 
eration defines a new, three-qubit operator [p3: 



\Ji\ 



Sj—= (x^ Si \Si,Sj,-SjW)- 

\l Ci \ 



\Ji\s 



I Si, Sj, SjW/ 



(2.35) 



The algorithm for simulating the closed chain "ID spin 
glass" can now be written as: 



rpj _ C\\Jn\ rpi 



(2.36) 



To prove that this formula indeed yields the correct ther- 
modynamic weights, we may employ the result for the 
open-chain case, which was proved in the previous sec- 
tion, and include an extra state for the workbit. In the 
calculation to follow, a ""' on the symbols indicates a 
sum over all spins except s±, and: 



JV 



<I> 



N 



(-irn 



(2.37) 



i=l 



Now: 



Tg\{-} N )\w = -i> = niSI^K-}*}!-) = 

2 N -1 

<iE ^^ [{sw i{s},)k=-i> = 



v=o 



— -^fx- H ^^\{s} y ){x^ s -\w = Sl )+x-^ s -\, 



"Si, 



v=o 



upon collecting terms with equal \w): 



<Z>n [^-^K-^'-'^W^-I^I^K-l, s 2 , s N } y ) + 



V CN Z N v 

(x- H N^ S ^--' S ^ X -^ SN \{l,S 2 ,...,S N }y) + 
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X[{-M2,..., S »} a ] J l*N|{_i iS2i ... i ^} !) )j |_) 



2™-l 



®nx~ h » [{s}v] \{ 8 }v) (x sjvsi|Jn| |+) + X - SNS1 



Jn\ 



cn a" n 



W E E ^^ HN[{s}sl ^ JNSNSl |{4y)|sign(Jiv)) 

^ AT . _ . 



sign(,/ N ) a=0 



2™-l 



sign(J N ) y=0 ^/ Z]^ 



(2.38) 



Since the amplitude-squared of the state |{s}y) is given 
by the thermodynamic weight of the corresponding ID, 
closed-chain spin glass system, we have proved that the 
algorithm which includes a workbit works for closed 
boundary conditions. 

How should one interpret the "|sign( Jn))" in the last 
line of Eq.(2.3S)? The workbit is in the excited- or 
ground state according to whether sign( Jn) is positive or 
negative, respectively. That is, the state of the workbit 
is determined by whether the interaction between spins 
si and sat is ferro- or antiferromagnetic. However, in the 
simulation of Ising models one is interested in a specific 
set of bonds, so it is necessary to be able to determine the 
last bond-sign. This is especially important for higher di- 
mensional models, where every plaquette corresponds to 
a closed ID chain. For this reason we will next present a 
detailed analysis of the complexity associated with gener- 
ating a single plaquette. Before doing so, we note that a 
useful corollary follows from the calculation above, given 
that there was no special importance to the indices of the 
spins between which the loop was closed: 

Corollary 2 Closing the bond between Si and Sj using 

always produces a superposition with half the states 
having probabilities equal to the thermodynamic weight of 
the Hamiltonian with a ferromagnetic Jij , and the other 
half with antiferromagnetic Jy . 

Now, the simplest way of determining the last bond sign 



detail. 



Jn 

N,l,w 



This 



is to measure it, following the application of f2 
irreversible, nonunitary operation collapses the superpo- 
sition of \w) while leaving the superposed state of the 
spin-qubits intact. It randomly chooses between a ferro- 
or antiferromagnetic bond connecting s± and sn, i.e., 
chooses between sign( Jn) — ±1. Measurements are tan- 
tamount to classical operations, so in light of the com- 
ment at the beginning of this section, a measurement as a 
means of choosing the last bond sign leaves the algorithm 
in the classical realm. Therefore it may not come as a 
surprise that measurements will prove to be ineffective as 
a means of obtaining the desired last bond sign. Instead, 
one may employ interference in order to "erase" one of 
the last bonds, and leave only the desired one. In the 
following two sections we will discuss both procedures in 



E. Measurements as a Means of Choosing the 
Desired Bond Sign 



A measurement of the workbit collapses its superposi- 
tion into either the ferro- or antiferromagnetic state. As 
we now show, which last bond sign has the higher prob- 
ability depends on whether the resulting closed chain is 
frustrated or not. From the result of Eq.(2.38) we have 
for the ± J spin glass: 



^ 1 naLt n \{-} n )\-) 




-/3<[M H ]/2| {s}y) | sign(Jjv)); (2.39) 



where J = (Ji, J2, Jn) and | = J. Let us now de- 
fine the "partial partition functions" Z N N , i.e., for 
a ferromagnetic bond between s\ and sn, and Z N for 
an antiferromagnetic bond. The relative weight of the 
ferro- and antiferromagnetic subspaces can then be ex- 
pressed as: r = Z~^/Z N . Without loss of generality, let 
us assume from now on that an antiferromagnetic last 
bond results in a frustrated system (i.e., the total num- 
ber of antiferromagnetic bonds is odd), and vice versa. 
Then r determines the relative probability of obtaining 
a frustrated or unfrustrated chain as a result of the mea- 
surement on the workbit. It is intuitively clear that the 
spin configurations of a frustrated system will generally 
have a higher energy than those of the corresponding 
unfrustrated system, and one would thus expect to find 
> Z N . In ID this statement can be made exact, as 
we now show. Separating the last bond one finds: 



N 



E 

y=0 



^-^(H :l N [{s}y}-J N s N s 1 ) 



(2.40) 



which can be split into two terms, corresponding to 
si = sn and si — — sn- 
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6 N 



e -p( Hl x[W'i=**]- j >') 



Mr 



{'}„ 



-Jn) 



(2.41) 



Defining H° N [{s}} = H N [{s} st=SN ] and H N [{s}} 
H N [{s} Sl =-s N }i we can write this as: 



N 



(x 2J »Z° N 



— 2J N 7I \ 
X N J ; 



{»} 



(2.42) 



Using this: 



x 2J Z N 



7+ 

7 - T -2J%Q 



„-2J Z l 



N 



x^zy 



The "constrained" Hamiltonians, where sn 
be written as: 



(2.43) 
±si can 



N-2 



^ JiSiS i+ i — Jn^iSn-isn 



H N 



— — Jn-isn~is\ 

— Hft-i + Jn-iSn-iSi. 



(2.44) 

This allows one to break up the constrained partition 



functions in a manner similar to Eq.(2.41 



= ^ e - j 3(gj r _ 1 [{8}]-Jjy-iJjf-i8i) 

W 



y^ e -/'(^- I [{»}.i=. W - 1 ]-J*r-i) 

{4 



{4 

/■^,2Jjv-i 7O 



_ 1 ]+-^7V-l) 



„— 2Jjv-i 7I 



Similarly, 

7I /"„.— 2JjV-l 7O 1 2J N -l 7I \ 

iv v ^Ar-i+a; ^n-i) 



(2.45) 



(2.46) 



One may continue to split up the Hamiltonians as in 
Eq.( p.44 ). The general pattern is seen to be (0 < n < 



7O 9 I 2J N -n-l 7O 

^N-n — * \ x ^N-n-X 

7I _ ( rr — 2JN-n-l 7O 



„ — 2Jw-n-l 7I 



2Jjv-ti-i 7I 



ZAr_ n _ x ). (2.47) 



Together with the obvious initial condition for a pair of 
spins, Z% = 2x 2Jl , Z\ = 2x~ 2Jl , this defines a recur- 
sion relation which can be solved to yield for Eq.( 2.42| ) 
(N > 3): 



)JV-2 



N-l 



Z 



N 



where: 



,4'E\ •/ 



N-l 



1 l\ zZ 



2 (J+(-l) k + N - 1 J N ) 



f N (k). (2.48) 



fc=0 



f N (k)= x^^ Ji i ; M0) = 1- (2.49) 

il<22<---<ifc=fc 

The last function generates all possible different sums of 
the Jj's. For example, for N = 4 we find: 



Z+ _ [1] + X 4J [x 4 ' h + X 4J * + X 4J3 ] + [x 4 ^ + J ^ + ^(•/i + ./a) + X 4(J 2 +J 3 )^ + X 4J ^4(./ 1 + ,/ 2+ J3)] 

Z7 ~ \x 4J ] + \x 4 ^ + x 4J 2 + a; 4J 3l + X 4J \x 4 ( J ^+ J ^ + z 4 ( J i+ J 3) + x 4 ( J 2+ J ^] + \ X 4 ^+ J ^+ J ^] 



(2.50) 



It can be checked that this agrees with the result 
obtained directly from the corresponding Hamiltonians, 
and indeed, Eq.(2.48) can be proved to hold by induction 
|35||. The important point about Eq.( [2.48 ) is the differ- 



ence between Z^ and Z N . As can be seen in Eq.( 2.5C| ), 



there is an alternating ratio of x between groups of sim- 
and Z N . "Similar terms" here refers to 



ilar terms in Z N 



terms with the same number of Ji 's, enclosed in brackets 
in Eq.(2.50). Since there is a one-to-one correspondence 
of this type, we may utilize the elementary inequality |pq| 



< 



J N 



Z 



, 4|J 

< x 11 = max 



(x 4J ,x- 4J ) 



(2.52) 



N 



The up per b ound is approached as T — > 0. To see this, 
use Eq. ( 2.48| ) to express the ratio of unfrustrated to frus- 
trated partition functions as: 



z + N _ Y. N k Jx 2J(1 - { - 1)k+N) fN{k) 



(2.53) 



mm — 



to obtain that: 



/ A I —4 / 

mm [x , x 



< 



) = 



= X -MJ\ 



< max ( — 

bi 



(2.51) 



Since Jn < results in a frustrated system, when N is 
odd there must be an even number k* (zero included) 
of ferro- and even number of antiferromagnetic bonds 
among the first N — 1 Jj's. When N is even, there must 
be an odd number k* of ferro- and even number of an- 
tiferromagnetic bonds among the first N—l Jj's. Now, 
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as T — > the do minan t term in both numerator and de- 
nominator of Eq.( 2.53| ) will be the one which has all the 
positive Ji's (such a term exists, since /jv(fc) generates 
all combinations of Ji's). When T is sufficiently low, all 
other terms become negligible, since they are smaller by 



at least 



„4J 



of Eq. ( |2.48 ) it suffices to consider that of the dominant 
terms. As can be seen from the expression for /jv(fc), k 
counts how many J's appear in every term. Thus for 
odd N the dominant term is generated when k = k* (is 
even), whereas for even N, the dominant term results 
when k = k* (is odd). In both cases k + N is odd. But 
this means that 2( J— (— l) k + N J N ) is zero in the frustrated 
( Jn = — J) case, 4J in the unfrustrated case. Thus the 
dominant term in the unfrustrated case is always greater 
by x iJ than that of the frust rated case. This proves that 
the upper bound in Eq.( 2.52 ) is indeed reached as T^O. 
The implication is that an algorithm which attempts to 
generate frustrated isolated plaquettes by measurement 
alone, will have to try an average of ~ x 4J times before 
succeeding, and this number grows exponentially as the 
temperature is lowered. 

One might be tempted to try to correct a "wrong" pla- 
quette instead of "discarding" it. However, it turns out 
that any correction procedure has probability less than 
1 of succeeding. For example, suppose one is interested 
in the frustrated case, but the measurement yielded an 
unfrustrated plaquette. The problem is then that 
includes the bond Jn with the wrong sign. One could 
imagine several strategies to "undo" this, which all start 
with a new workbit w' . For example, one could employ 

a "biased random walk" procedure: one redoes fij^iL'i 
measures again, etc. The hope is that the resulting sum 
of Jjv's with random signs will at one point add up to 
produce the originally desired ferromagnetic bond. The 
probability for this to happen is equal to the probabil- 
ity of return of some biased random walk where the bias 
increases with the distance from the origin. For the un- 
biased random walker in ID and 2D this probability is 
1, but the waiting time is infinite J57|| . For the biased 
random walker the probability of return turns out to be 
zero. This means that a "random walker" correction pro- 
cedure cannot guarantee the desired result, at constant 
temperature. We are not aware of any other, more suc- 
cessful procedure. What about turning a frustrated last 
bond into an unfrustrated one? Even this has probabil- 
ity less than 1, as a consideration of a three-spin system 
will illustrate: Suppose one chooses J± = J2 = J and af- 
ter measuring w one finds J3 = — J. One might then 
hope to correct this frustrated system by redoing fig J ' / . 
A low T analysis will suffice to show that this will not 
correct the error. At low T the dominant spin config- 
uration will be that with the lowest energy E m i n . It 
is easily checked that both w' = ±1 yield E m i n = — 2J 
and thus have equal probability. The case w' = 1 cor- 
responds to a ferromagnetic last bond and therefore cor- 
rects the Hamiltonian. However, in the equally probable 



opposite case the Hamiltonian now includes a last bond 
of strength J3 = —2 J. If this had been the result of the 
correction procedure, one would be facing a (wrongly) bi- 
ased random walk again, since a new measurement would 
either result in J3 = J or J3 = —3 J, with corresponding 



Thus to understand the low T behavior E n 



-J and E„ 



-3 J. The latter is more probable 



by a factor e 2/3J , so the correction fails. 

These arguments show that procedures using only com- 
binations of superposition and measurement, have no 
control over the type of plaquette they generate. In the 
next section we will introduce a procedure which does 
have this feature. 



F. Using Interference to Close an Isolated Plaquette 

After closing the last bond, the quantum register is in 
a superposition corresponding to a frustr ated (w = — 1) 
and unfrustrated (w — 1) plaquette [Eq.( 2.38| )] (assum- 
ing, without loss of generality, that there is an even num- 
ber of antifcrromagnetic bonds among the first N — 1 
bonds): 



IV') = T*\{-} N )\w = -I) = |V;-> + |V+>, 



(2.54) 



where 



IV'd 



<S>NX- H ±^\±)\{s} y ), (2.55) 



and where H± corresponds to Jjv = ±J. We have inten- 
tionally written the workbit first, so that in the binary 
representation of the spin configurations by the quantum 
register, the first 2 N states correspond to the frustrated 
configurations (w = —1) and the last 2^ correspond to 
the unfrustrated configurations (tu = +l). As was shown 
in the previous section, achieving control over the plaque- 
tte type cannot be done by measurement alone. However, 
one may try to employ interference in order to "erase" 
one of the subspaces, thus leaving only the desired pla- 
quette type. To see this, consider the wavefunction \ip) of 
the quantum register as a vector of length 2^ N+1 \ with 
the first 2 N entries corresponding to \tp~), and the last 
2 N entries corresponding to \i[>+). Within each such sub- 
space, the entries run over all possible spin configurations 
y = 0...2 N - 1. Clearly =1 by unitarity of T$ . We 

now seek a new unitary transformation U such that: 



U±\il>) = 



1 



± ■ 



(2.56) 



Thus U rotates the superposed quantum register state 
into a state representing either the frustrated or the 
unfrustrated configurations. That U exists is clear since 
it takes one norm 1 vector into another. Furthermore it 
clearly mixes different spin configurations, i.e., creates 
interferences. The problem is to find U , given that it is a 
2(iv+i) , 2i N + 1 ) matrix of coefficients that depend on J . 
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In other words, one needs to know the Gibbs distribution 
of the plaquette as input in order to find U\ This might 
seem to defeat the purpose of the algorithm altogether, 
but not so in view of the need to construct plaquettes 
with given bonds for the 2D and 3D Ising problems. In 
these higher dimensional cases it is very useful to know 
how to produce small plaquettes of, e.g., 3 or 4 spins, 
for the triangular and square lattices respectively. Thus 
in these cases, or indeed for any plaquette of N spins, 
one could calculate in advance the Gibbs distribution, 
find U, and use it in the construction of a lattice. We 
will deal with the quest ion of "integration" of a plaque- 
tte into a lattice in Sec. Ill B. Here we give the general 



n-1 



2n-l 



(classical) algorithm for the construction of U for any N, 
and explicitly solve the cases N = 3,4. 



1. General Construction of N -Spin Interference Operator 

Let {ei}?" x , n = 2 N , be the standard basis of vectors 
for TZ 2n , with a 1 at position i, zeroes elsewhere. Con- 
sider a normalized real vector v of length 2n (representing 
|V>)) and another vector w composed of v's upper half, 
also normalized: w = (i>i, t>2> •••) v n , 0, 0, 0)/a_ where 

ct_ = (X^r=i^) 1 ^ 2 - Here w corresponds to \ip-). We 
seek a construction by two-qubit operations of a unitary 
matrix U_ such that U_v = w [as in Eq.(2.56)]. The so- 
lution to this problem is within the theory of generators 
of SO(n), as outlined in Ref. j580, in terms of generalized 
Euler angles. Since in our case v and w are real vec- 
tors, it suffices to find an orthogonal U_. The idea is to 
first rotate w so that it coincides with e„ , and then solve 
the easier problem of finding the transformed U_ which 
rotates the transformed v into e n . More explicitly, sup- 



pose one has found an orthogonal matrix G_ satisfying: 



w = G_ e„. Clearly G^' is composed of two blocks, an 
upper block G'^ and a lower one 7„, the n-n identity 

(2) 

matrix. The transformed equation is then: G_ v = e„, 
with GL 2) = {G[ 



v' = (G W )- 1 v=(G (1) )- 1 («-W+ >T v t ei) 



i—n-\-l 



(0,0,..,a-,V n+ i,...,V2n)- 



(2.57) 



Since only the last n + l coordinates of v' are non-zero, 

(2) 

G_ is composed of an upper block 7 n _i and a lower 

block G'i 2) which we need to find along with G_ . Hav- 
ing found these, the solution can be written as: 

U_ = G [ 1 ) G^ ) {G W )- 1 . (2.58) 
Following Ref. |58|], let us write: 



c a) = n 9i@i), gl 2) = n gM 



(2.59) 



i=i 



where gi(9i), is a rotation by 6i in the plane spanned by 
{e t ,e l+1 ) ofTl 2n : 



( h-\ 



9i{0i) 



cos t>i sin t 
- sin 0i cos ( 



(2.60) 



V 



Application of gL 1 ' to e„ results in the following set of 
equations: 

W n = COS^-i 

w n -i = sin n —x cos B n —i 

w n -2 = sin^j-i sin6'„_2 cos#„_ 3 



W2 = sin u n -\ sin W„_2 • • ■ sin c/2 cos V\ 
w\ = sin0„_i sin n _2 1 ■ 1 sin 62 sin^i, 

with solution (k = 1, n — 1): 

Vk+l 



(2.61) 



COSPfc = 

sin. 6k = 



Tk+l 



1/2 



(2.62) 



Similarly, application of (G^) 1 to e n results in 



V n = COS ti n 

v' n+1 = sin0„ cos6»„ + i 

v' n+2 = sin6 n sin0 n+ i cos0, 1+ 2 

V 2n-1 = Sm0„ Sin6» n+ i • • • Sm0 2 n-2 COS02ri-l 

v' 2n = sm9 n sin6>„+i • • • sin6»2„-2 sin 9 2n -i, (2.63) 
yielding (k = n, 2n — 1): 



cos Ok = — - 

Tk 

■ a r k+i . a r k+ i 



2n 



1/2 



(2.64) 



The case w = (0, 0, v n+ i, v 2 n)/c+, with a+ = 
1/2 

(X)i=i v i) > corresponds to 1^+), and we need to find 
an orthogonal U+ such that C/+v = w. It can be 
seen by repeating the arguments above that we now 
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require w = GT'e n -|_x, whence G^P has /„ as its up- 
per block, and v' = G^v = (v\, ...v n ,a+,0, ...,0) sat- 
isfies the transformed equation v' = G^ e„+i, where 

G^ 2) = G^V+^G^V 1 has J n _i as its lower block. 
Writing accordingly 



7V-1 



G 



(i) 



n am 



G 



(2) 



(2.65) 



i— n+l 



leads to equations very similar to Eqs.(2.61-2.64): one 
needs to replace Wi and Vi by everywhere in Eqs.(2.61 



2.62| ), as well as allow k to range from 1 to n. In Eq.(2.63) 
one should replace v[ by Wi , whereas in Eq. ( 2.64 ) k should 



range from n + 1 to 2n — 1 and v[ replaced by Vi . With 
these replacements it follows that the 9i are identical for 
U- and U+ , except for 9 n , for which sin = cos 9~ . The 
interference matrix is given by: 



(G^y'G^G^ 



(2.66) 



This, together with Eqs.(2.57),(2.58), uniquely solves our 
problem. It remains to be shown explicitly how U± can 
be written in terms of one and two-qubit operations on 
the original spin-qubit register. 

We note that the gi(9i) are identity matrices except 
for 2 • 2 blocks. Consider the representation of the ba- 
sis {e^} by s' = {\s[, s' N )} y ,S~ where the successive 
register states differ by a single qubit flip. In this basis 
application of each gi(8i) does not entangle states differ- 
ing by more than one qubit, so the corresponding gi(9i) 
is a single-qubit operator: 



i (e i )|a i = -l) = cos0 i |-)+sin0 i |+) 
i(0i)\s t =+l) = -wx9i\-) + cos 6i\-\ 



(2.67) 



The remaining problem is thus seen to be the transforma- 
tion from the original "binary" basis {|si, sn)}^ = q 1 
to s' and back. This is easily accomplished by two-qubit 
operations using the well-known classical Gray code |^9| . 
The Gray-to-binary transformation is accomplished by 
successive XORs starting from the last qubi t [these dif- 
fer somewhat from the definition in Eq.( 2.33| ), hence the 
different notation]: 



X 'i,j\Si 



(2.68) 



where an extra workbit w x represent sq = — 1. For 
example, for the two successive binary s basis states 
\w x , St, s 2 , s 3 , s 4 ) = ]-,-, +,+,+) and |-, +,-,-,-) we 
find after application of Yl i=0 X- | — , — , +, — , — ) and 
| — ,+,+,—, — ) respectively, which indeed differ by only 
a single qubit. Furthermore, clearly X'^- —X- j so the 
binary-to-Gray transformation is accomplished by run- 
ning the same sequence of XORs in reverse order (start- 
ing from the extra workbit). We can now finally write 
down the full interference transformation. Let: 



X N — X W:i; W X w l x^ i+1 , 



(2.69) 



and: 



2"-l 



i=l 



i=2 N 



(2.70) 



=2 N + 1 



i=l 



The inverse operators are obtained by reversing the order 
of the products and negating all angles. Then: 

\w x )\tP.) = X'^G^G^iG^r'X'^M 

\w x )ty + ) =X^\G^Y\G ( + ) r 1 G ( + ) X' N \w x M (2.71) 

expresses the interference transformation leaving only 
the register states corresponding to frustrated or unfrus- 
trated spin configurations. 



2. Solution for N = 3, 4 Spin System 

We now employ the above formalism in order to ex- 
plicitly solve for the interference transformations of the 
three- and four-spin systems, corresponding to the ele- 
mentary cells of the triangular and square lattices respec- 
tively. First it should be noted that for the ± J model in 
ID, for a closed Ising chain of given length, the spectra 
of all frustrated realizations of bond choices are identical, 
and so are the spectra of all unfrustrated realizations. To 
see this, consider a specific spin configuration and real- 
ization of bonds with energy E, and suppose one changes 
the sign of some arbitrary pair of bonds J rn , J n , m < n. 
This operation does not change the frustration of the 
chain, since this is determined by the parity of antifer- 
romagnetic bonds. But by flipping all spins s m +i, ...,s„ 
once again the energy is E, since the flipping of s m +i 
and s n undoes the change in sign of J m and J n respec- 
tively, and all other spin flips occur in pairs that share 
a bond and thus cancel. So for every spin configuration 
and choice of bonds, there is another spin configuration 
with the same energy in a realization with the same frus- 
tration but different bonds. Clearly the mapping above 
is one-to-one, so that indeed the spectrum is identical 
for all bond realizations with the same frustration. Re- 
turning to the N — 3,4 spin systems, we are at liberty 
to consider, e.g., the case where all bonds but the last 
are ferromagnetic. The superposition of J/v into a ferro- 
and antiferromagnetic bond will represent all other un- 
frustrated and fr ustra te d bo nd realizations, respectively. 
Solution of Eqs. fl2.62 ), (2.64) yields the following result 
for the transformation from the superposition to the frus- 
trated or unfrustrated subspaces. Let: 



14 



D.A. Lidar & O. Biham 



fk,l,m — 



„4Jfc 



VTT" 



r 8.7 



hk,l,m,n — 



\Jl + mx 8J + nx leJ 



(2.72) 



Then, writing C = cos, the angles can be expressed as: 



N = 3: 
COi = 

ce 15 = 

ce- = 



i 

V2 



1 

V2 



ce 14 = ^ 

(/0, 1 .3(l+^ 4J ) 3/2 )' 



C#3 = — /o,l,3 
C013 = /l,l,3 

= sin #t 



C&l = /l,l,4 6*05 = /l,l,5 C*^ = — /l,l,6 C*#7 = — /o,2,l 

C012 = — /(),4,1 C#n = — /o,5,l C$10 = /o,6,l C#9 = /l,6,2 



AT = 4: 
C0i = 

ce 3 = 
C0 S = 
ce 8 = 

COu = 

C013 = 

sin 9± 6 
C9 17 = 

C6>19 = 
C025 = 
C026 = 



73 CT2 = -71 

— fro, 1,3,0 C^4 = frl, 1,4,0 
^0, 2, 4,0 C8e = —fro, 3,4,0 

— frl, 4, 5,0 C6g = — frl,4,6,0 
: fro, 5, 7,0 C$12 = — frl,5,8,0 
: — fr(), 6, 8,0 C$14 = fr-0, 7, 8,0 

(fro,l,6,l(l+X 4J ) 2 ) _1 = C0+ 

: — fr-2,2,12,2 C$18 = — frl,2,12,l 

: frl,2,ll,l C02O = frl,2,10,l C#21 = 
: fr-2, 1,6,1 

: hi 1,6,0 C$27 = —frl, 1,5,0 C#28 = 



C87 = —fr-0, 4, 4,0 
C#10 = frl, 4, 7,0 

C#15 = fro, 8, 8,0 



-frl, 2, 9,1 C$22 = —frl, 2, 8,1 C#23 = frl, 2, 7,1 C$24 = fr-0, 2, 6,1 
"frl, 1,4,0 C$29 = frl, 1,3,0 C830 = frl, 1,2,0 C831 = —frl, 1,1,0 
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The regularity and similarity between terms in the same 
column (cos Oi , cos 9j with i+j — 16) for N = 3 is notewor- 
thy. For the four-spin system there is a similarity between 
terms in the same row, but we find a less regular solution. 
Let us remind the reader at this point the motivation for 
intr oduc ing the above transformations. We showed in 
Sec. HE that the average number of attempts needed to 
generate a plaquette of given type using only superposi- 
tions and no interference, grows as x 4J = exp(2 J/Zc^T) 
with the temperature. Using the interference transfor- 
mations, the cost is 0(1) in the plaquette size, and inde- 
pendent of the temperature. We are now finally ready to 
discuss more interesting Ising models, in 2D and above. 



III. HIGHER DIMENSIONAL ISING MODELS 

The "ID Ising spin glass" is rather trivial, and the 
more interesting models are the higher-dimensional ones, 
where connectivity plays an important role. As an intro- 
duction to the schemes we will need to employ in dealing 
with the 2D and 3D cases, consider first the "infinite- 
dimensional" Bcthc lattice, which has no closed loops. 



A. The Bethe Lattice Case 

Consider a binary Bethe lattice, i.e., with spins located 
at the vertices of a binary tree (Fig.[y). The spin glass 



Hamiltonian for a if -level deep tree can be written as: 

K 2 k 

h k = ~X]X! G ( fe -UV2l)(M) s (fc-UV2"l) s (M)- C 3 - 1 ) 



According to the general recipe of Sec. II A, the quantum 
operator for calculating the weights of configurations in 
this system is (for simplicity of notation we shall suppress 
the indices on G where they are already indicated by S) : 



If 



nn^ 

h=Ki=\ 



rv2i),(fc,i) 



R 



(0,1)- 



(3-2) 



In order to prove that Tk correctly calculates the proba- 
bilities of the spin glass Ising model on the Bethe tree, we 
need to show (I) that it does not matter in which order 
we connect the spins occupying vertices one level deeper 
then their common originator, and (II) that a ID chain 
which splits at its end into two branches is correctly de- 
scribed. For, (I) allows us to perform the first branching 
in the tree [from spin (0, 1)], and (II) [combined with (I)] 
allows us to build up the tree recursively from any ex- 
isting end point. In particular, the order described in 
Eq. ( |3.2| ) will be valid. Starting with (I) , we need to show 
that [S G , Sf^\ — (the indices i,j,k are shorthand for 
the double indices employed above). Now, 



qGqG 

on Oik 



1 



Si , S 



Cik 

1 



Sij [: 



s k x G ^ s > 

1 

s k x {G ' k - 



3,sk) 

-G thSi 



Si , Sj, S k ) 



-GijSi 



S k X 



G ikSi 



S{ , S.4 



'l! °3 1 



8 k) 



r G ijSi 



J ' 



Sk) 



Sj,Sk) +SjX 



3 i) °j > 



GijSi 



-(G <fc +Gy) Sj | 



,Sk) 



.(G„ 



-Sk) + SkSjX 



(G ik +Gi 



1 1 Si , Sj , ■ 



, Sfc> + 



Sk) 



(3.3) 



But on the other hand, by exchanging j and k everywhere in the last line, we obtain: 

&ik $ij 1 3 i ; s j i s k) = 

x -(G l]+ G tk )s t | Si) Sfe] + SkX (Gik-G v ) Si ^ _ Sfc; + 



yfcihi 



^- G ^Si,S k ,- Sj )+S j S k X^ +G ^\ f 



(3.4) 



The order of the qubits in the kets is immaterial, so that 
by comparing the two results we find that indeed 



[S^, S iki 



0. 



This 
Next 



is 

we 



indicated graphically 
prove (II) above, 



(3.5) 

in Fig.|(a). 
namely that 



s n,n+2 S n,n+i t n\{-}n)\-)n+i\-)n+2 (where due to 



Eq. ( |3.5| ) we may exchange the order of S G N+2 and 



S 



G 

JV.JV+l/ 



yields the correct thermodynamic weight for 



the Hamiltonian 



1G 
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N-l 



H — ^ Gi^ + iSiS i+ i + Gn^n + is n s n+ i 
i=i 

Sn,N+2^N,N+1 T N |{-}jv)|-)iV+l|-)jV+2 
2 N -1 

^ n ^N,N+2^N,N+i\Vn)\~ ) ~ ) = 



+ Gn.N+2SnSn+2- 



(3.6) 



Using the results of Eqs.( 2.25 ), (3.3) we have: 



v=o 



2"-l 

\ ^ j^G I— \ l x ~s N (G NtN + i+G N ,N + 2) I \ j,sn(Gn,n + i — Gn,N-\-2) \ 



y/CN,N+2CN,N+i. ir (1 

j- s Af (GtV,7V + 2 — GjV,JV + l) J _j_ j> s JV (GjV,iV + l+GjV,iV + 2) ~t~)J 

1 



y/CN,N+2CN,N+l ^ 



SAT+lSAr+22; 



G J v,iv + i3«3iv + i+G N , N + 2SN;iJ v + 2|y Ar ^| SAr+ljSjv+ ^ = 



^CAr :A , +2 CAr iA r + l ^ AjjG 



1 1 

c -P(H^\is}v\—GN + lSNS N +l—GN + 2SNSN + 2)^y^^\ 



(3.7) 



which is the desired result. It should be noted that since 
[S^S^] — 0, any number of S operators with a com- 
mon starting point i will commute in pairs. Therefore 
there is nothing special about the binary Bethe tree, and 
we can equally well app ly our algorithm, after a proper 
modification of Eq.( |3.2| ), to a Bethe tree with any kind 
of branching. 



B. 2D Ising Model 

As was demonstrated in the ID case, the key to being 
able to close loops is the creation of a superposition in 
frond-space by using a workbit whose state is compared 
with the spin with which the loop is closed. Choosing 
a specific bond-sign is then accomplished by an inter- 
ference transformation which eliminates one of the bond 
subspaces. We now extend these ideas in order to present 
an algorithm for simulating 2D Ising spin systems. Ide- 
ally one would like to have an algorithm which can ex- 
actly calculate the thermodynamic weights of an arbi- 
trary given spin glass Hamiltonian: 



H 



Jij s i s j 

{id) 



J. 



(3.8) 



However, since the interference transformations intro- 
duced in the ID case require as input the thermodynamic 
weights, we cannot hope to deal with an arbitrary Hamil- 
tonian. Instead, as will be shown here, the class of spin 
glasses which can be dealt with by the present algorithm 
is that with predetermined plaquettes of finite size. In 
other words, by using interfence transformations one can 
construct isolated plaquettes of any desired (finite) size 
and composition (of bonds), and these plaquettes can 



then be connected together. This creates new plaquettes, 
with random bond signs. Thus the algorithm cannot pro- 
vide complete control over the bond composition of the 
Ising model it is used to simulate, but the resulting class 
of partially random-bond systems is huge (exponentially 
large in the number of bonds). Furthermore, by lowering 
the temperature one can increase the probability of gen- 
erating only unfrustrated plaquettes connecting the pre- 
fabricated ones. For example, choosing prefabricated un- 
frustrated plaquettes will generate low-temperature sim- 
ulations of unfrustrated Ising models with very few de- 
fects. That is, if Nd denotes the number of defect (i.e., 
frustrated) plaquettes, then: 



Nd 
N 



„-4J 



(3.9) 



We turn next to demonstrating how isolated plaquettes 
can be connected together to form a lattice. 



1. Allowed Algorithms 

"Hooking up" isolated plaquettes will require connect- 
ing spins by operators, all of which will eventually 
have to share lattice points in pairs (or more), as shown 
in Fig||. Corollary || assures that bonds can always be 
closed using Q, so as to produce the correct superposi- 
tion. Since the order in which the lattice is constructed 
might appear to be important, the question of commu- 
tation of the various operators naturally arises. In this 
section we will deal with this in some detail. As for pairs 
of fl operators, all possible combinations commute [see 
Fig.|b(I)-(IV)]: 



\J\ 1 = 

ikw2 -I 
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[^lyLi ' ^LjL 2 ] = ^ ' [^Ijli ' ^wlu 2 ] = ^- (3.10) first of these relations (we drop the normalization factors 

, ,, . , ,. . , ,, and set J = 1 for notational simplicity): 

Wc demonstrate the calculation required to prove the 

^jkw 2 ^ijw 1 \Si,Sj,Sk,Wl,W2) = 

tt 3 kw 2 [~Sj (x Sz \s i ,S j ,S k ,-S :j Wi,W2} - WiX~ S *\Si, Sj, S k , SjWi,W 2 ))] = 

S k Sj [x Sj+Sl \s i ,S :j ,S k ,-S : )W 1 ,-SkW2) -W2X~ S: l +Si \Si,Sj,Sk,-SjWl,SkW2) + 

-WiX S3 ~ Si \Si,, Sj,S k , SjWi,-S k W2) + WlW2^ SrS " | Sj, Sj, S k , SjWl, s k w 2 }] , 

whereas on the other hand: 

^■ijw 1 ^jkw 2 \Si,Sj,S k ,Wi,W2) = 

&ij Wl [—s k (x Sj \si,Sj, Sk, wi, -s k w 2 ) - W2X~ S] \si, Sj, s k ,wi,s k w 2 ))] = 
SjS k [x Si+Sj \s i ,s : j,s k ,-s : jW 1 ,-s k W2} - wix~ Si+s i\si,s 3 ,s k ,s j w 1 ,-s k w 2 ) + 

-W 2 X S *~ Si \Si, Sj,S k , SjWl, S k W 2 ) + WiW 2 X~ Si ' Sj \Si,Sj, S k , SjWl, S k W 2 }] ■ 

It can be verified that the last lines in these two calculations arc identical, proving the first commutation relation. 
Next we consider combinations of <S* and f2. They all commute except one: 

[s44L] = o; [^Lt] = o; [<^,nif]j = o; [s^n^j, (3.11) 

[S^n[%}\ SuSj , Sk ,w) = ~^L=x Js > (x Js > ■ wx Js > ) \ Si ,- Sj ,Sk) [\s jW ) - \ - Sj w)} (3.12) 



/CiCfc 

We demonstrate this, again taking J = 1 and dropping normalization: 

^kjw Sij | Si , Sj , Sk, W) 

tt k j W [x~ Si \si,Sj,Sk,w) + SjX Si \si,-Sj,Sk,w)] = 

-Sj (x Sk ~ Si \si,Sj,Sk,-Sjw) - wx~ Sk ~ Si \si, Sj,s k ,s k w) + SjX Sk+Si \si, -Sj, s k , Sjiv) - ws j x~ Sk+Si \s l ,-s J ,s k ,~s j w}) 
whereas: 

Sij Qfsjw | $i , Sj , S k , 

-SjSij [x Sk \Si, Sj, Sk, -SjW) - WX~ Sk \Si, Sj, Sk, SjW)] = 



- Si \x- s * +Sk 



\Si,Sj,Sk,-SjW) + SjX S ' +Sk \si,-Sj, S k ,-SjW) -WX S * Sk \Si, Sj, S k , SjW) - WS X S% Sk \Si,-Sj,Sk,SjW}] , 



which is different from the previous result in the sec- 
ond and fourth terms. The reason that the relation in 
Eq.(3.12) does not commute is that only ilkjwSij pro- 



duces the right result. In the opposite order the same 
problem arises as when one tries to close a loop with S"s 
only. This relation is depicted graphically in Fig.H( d). 

Given the commutation relations in Eqs.( |3.10| ) and 
(3.11) one is essentially free to connect the isolated pla- 
quettes in any order, as we prove next. 



2. Proof of the Algorithm 

For simplicity we consider the case of a square lat- 
tice, where one has prepared a set of 2x2 plaquettes and 
placed them with equal spacing on an N by M lattice 
(Fig||). Denote the Hamiltonian for this system by Hq . 
Clearly the maximum density of non-overlapping prefab- 



ricated 2x2 plaquettes which can be achieved is 1 /4. This 
can be increased by using larger plaquettes, at the price 
of increasing complexity in their fabrication. The prob- 
lem is now to connect plaquettes, and from Corollary |^ 
this can be done with f2 operators which connect two oc- 
cupied lattice points. The geometries which may arise in 
connecting plaquettes are summarized in FigJ|(b-d) . The 
commutation relations of the previous section show that 
the only problem can arise in the geometry depicted in 
Fig.||(d). However, as long as S is applied before f2, the 
outcome is a correct superposition provided k Figj^(d) is 
the index of an occupied site. The commutation relations 
assure that in any order in which the O's are applied the 
lattice is generated correctly. We may thus assume that 
some arbitrary sequence of fi's has been applied. We as- 
sume further that workbits corresponding to new bonds 
are measured after application of fi, so that they are no 
longer in a superposition and a bond with random sign 
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has become integrated into the lattice. It will then suffice 
to prove that introducing a bond at an arbitrary location 
in the existing lattice produces the correct superposition. 
Let us proceed by induction, and assume that some par- 
tial set of all bonds has been closed by the algorithm. 
These bonds can be either horizontal or vertical, and for 
dcfinitcncss we will assume that they are always closed 
rightwards or upwards. Let us denote the set of K verti- 
cal bonds by {(n,m)}K and the K' horizontal bonds by 
{(n', m')}K'- The Hamiltonian for this set is: 



It should be noted that every intermediate stage in the 
construction of the full 2D lattice corresponds to a diluted 
2D lattice, with bonds specified deterministically by , 
and other bonds being either ferro- or antiferromagnetic 
depending on the result of measuring the workbits. Thus 
it follows that the algorithm can be used to simulate di- 
luted Ising models as well, with the same restrictions 
applying to the generality of the class of these models 
as specified above, i.e., the size of frustrated plaquettes 
must be finite. 



u 3 



HZ 



^ ^ J (n ) m) 4 'H,JTi 5 n,m+l 
{(n,m)}if 



E 

{(n / ,m , )} K / 



(3.13) 



with a corresponding operator: 



T (K,K>) — IT { M ™n,™^ l [ n [m)(l,m+l)w„ 



{(n,m)} K 



\ J 't I M 

' (n' ,m')(n' -\-l,m')w n 



T, 



(3.14) 



represents the measurement of workbit 



n 

{{n',m')} K i 
Here M Wn< 

). The notation Y[{( n . m )} K i m P nes a product over 
all members of the set {(n, m)}jf, in any order, and sim- 
ilarly for the horizontal bonds (guaranteed correct due 
to the commutation relations). Since the measurements 
have taken place, the workbits are no longer in a super- 
position after T? K K ,-. has been applied. Therefore the 
induction hypothesis is: 

T (K,K>)\{-} N m)\{-} K+K') = 
1 



2" " -1 

E 

y=o 



Z(K,K') 



IW*> 



(3.15) 



To see this it is helpful to consider for a moment the 
situation that would arise without the intermediate mea- 
surements: for every new (horizontal or vertical) bond 
|^n,m| one closes, one needs to introduce a new workbit 
w n ^ rn . A total of K + K 1 bonds thus requires a workbit- 
register |{— }k+k') as above. After \ J n , m \ is closed, the 
workbit is in a superposition of states corresponding to 
sign(J„. m ) = ±1. This superposition is destroyed after 
the measurements, leaving the spi n-qub its' superposition 
intact. This is the content of Eq.( |3.15[ ). 

The induction proof now requires showing that clos- 
ing an additional bond, say |J P , g |, produces the correct 
superposition jSOt l - Th e calculation is essentially identi- 
cal to that of Eq.( 2.38| ). It is easily checked that (i) this 



calculation is insensitive to whether | J Pi9 1 corresponds to 
a horizontal or vertical bond, and (ii) the appearance of 
a double index for every spin does not make any differ- 
ence. There is thus no need to repeat the calculation 
of Eq.(2.38), and we conclude the 2D algorithm to be 
proved. 



3. Control of Bond Sign 

One may wonder whether the lack of control over the 
bond sign is special to the essentially ID situation of 
closing an isolated plaquette. In fact the same prob- 
lem arises whenever a bond is closed using f2, as we 
argue next. Let us suppose that the algorithm has 
correctly produced the thermodynamic weights of the 
Hamiltonian Hq = —"Y^uj\JijSi8j, describing part of 
the full 2D problem, and excluding in particular the 
bond J n m- Associated with Hq is a partition function 
z o = J2{s} ex P(-P H o[{s}})- Let p nm = Pr[sign(J nm ) = 
1) be the probability of a ferromagnetic bond, q nm = 
Pr[sign(J„ m ) = —1] that of an antiferromagnetic bond. 
When the new bond J nm is included, we have for the 
ratio of these probabilities: 



Pnm 
Qnm 



a}' 



-l3H Q [{s}} 0\J nm \ SnSr) 



(3.16) 



This expression can be bounded from above and below 
by replacing s„s m in exp(±/3| J m „|s„s m ) by +1 or -1: 



|J " ml E' 



-PH„[{s}] 



< 



^ e -l3H a [{s}} e -l3\J nm \s n s m ^^ e ~(3H Q [{ s }} e l3\J nm \ 8n s n 

M {s} 



-0H o [{s}] 



Inserting this into Eq.(3.16) yields: 



-41 J n 



Q_nrn 



(3.17) 



(3.18) 



This resu lt is very similar to that obtained in the ID 
case, Eq.( [2.52" ), and although we cannot perform an ex- 
plicit calculation here, the implications are likely to be 
the same. Namely, that the upper bound is approached 
as T— > so that a frustrated plaquette becomes exp(2/3 J) 
less likely than an unfrustrated one. The main difference 
compared to the ID case is that now closing a single 
bond J nm may correspond to closing several plaquettes 
at once, which can only amplify the effect. Thus the is- 
sue of control of the bond sign is even more problematic 
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in 2D and above. Nevertheless, since isolated plaquettes 
can be "prefabricated" , the algorithm covers an exponen- 
tially large class of Ising models, as argued above. 



C. 3D Ising Model 

Three-dimensional Ising models are notoriously diffi- 
cult, both analytically and computationally. In particu- 
lar, whereas a large number of 2D models have yielded to 
analyses |3l| . only very few 3D systems have been solved 
analytically p3-p5[. Even so, all these models suffer from 
various shortcomings, such as having negative Boltzmann 
weights |32],|3| or being essentially 2D |64||35| . Computa- 
tionally, the plagues of 2D are of course multiplied in 3D, 
making our knowledge of 3D systems limited. Moreover, 
there are good reasons to suspect that calculations in 3D 
are fundamentally harder than in 2D. For example, some 
3D problems, such as finding the spin glass ground state, 
arc known to be NP-hard This is related to 

the existence of a freezing transition at T c > in 3D, in 
contrast to 2D, where T c — |6(|. Physically, this means 
that for T < T c the system can get stuck in a local min- 
imum in 3D (ergodicity breaking) and never reach the 
ground state. Thus it is extremely important to develop 
efficient algorithms for 3D systems. Let us now describe 
the extension to 3D of the 2D algorithm presented in the 
previous section. 



the presence of a magnetic field. The idea is to generalize 
the basic R and S operators introduced in the ID case. 
Let: 



1 



0,A. 



(a;** | ) SjiE - s *|+)) , (4.1) 



and 

S J ^\ Si , Sj ) = 

L= (x-V'''+W\8 i ,8 j ) + 8 j x J <'< +A '\a i ,-8 j )) , (4.2) 



c. . 
where 



c JiAi = 2cosh[/3(J, + S;Ay)]. 



(4.3) 



It is easily checked that i? Ai and S Ji,Aj are unitary. The 
loop-closing operator SI can be generalized accordingly: 

il \si,Sj,w) = X j!W S\^j' lAil \si,Sj,w) = 



l^l,|A 3 



(xWsi+\±A\ ShSj ,- SjW y 



wx-Q J *\« + W)\a i ,a j ,8 j w) 



A. Open Chain Case 



(4.4) 



1. Algorithm for 3D and Higher Dimensions 

The algorithm for the 3D case is a natural extension 
of that for 2D and presents no essentially new problems. 
One can either prepare 2D plaquettes and stack planes 
connected by S7 operat ors, o r prefabricate 3D cubes (us- 
ing the methods of Sec. II F) and connect those with Si's. 
In both cases Corollary |2| and the commutation proper- 
ties of the S7 operators guarantee that the correct super- 
position is produced, and the order in which plaquettes 
or cubes are hooked up does not matter. The only new 
feature is that more than two bonds can now emanate 
from the same site (for a cubic lattice). However, these 
present no problem since the algorithm is invariant to 
the order of creation of such bonds, as they all commute 
in pairs. In other words, the induction proof for the 2D 
case holds here as well, and the 3D algorithm is proved. 
Clearly, also the complexity remains the same, namely 
O(N) as long as random bonds between plaquettes are 
allowed. Moreover, the same argument holds for any di- 
mensional Ising system. 



IV. INCLUDING A MAGNETIC FIELD 

In this final section we show how the present algorithm 
can be extended in order to deal with the Ising model in 



In order to introduce an arbitrary magnetic field on 
every spin in an open chain geometry, consider the effect 
of applying R and S on a two-qubit register: 



n A. ' 



=a; JlslS2+AlS1+A2S2 |si,s 2 ), (4.5) 



0,Ai ' / J!,A 2 

Cj_ «1,»2 \ C s 



omitting some intermediate lines of by now familiar al- 
gebra. The coupling of Aj to the Sj in the exponent is 
like that of a magnetic field. However, the normalization 
factor also depends on s\, so it needs to be considered as 
well: 

J, A _ 



log c: 

i [log(> J + A )+e^ J + A >) + 
log (V 

£ [log (e^+ A ) 
log (e^ J ~ 



/3(J-A) + e -/3(J-A) 



-/3(J+A) 



-A) , e -0(J-A) 



- log (4cosh[/3(J + A)] cosh[/3(J - A)]) + 

a / cosh[/3(J + A)] 
2 ° g U°sh[/?( J - A)] 
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so that: 



B. Closed Chain Case 



exp | -- 10- r 



J, A 



J. A J. A 

C_ c: 



— 1/4 1 , I J, A / J, AN 

x -2fjS\og( C+ /c_ 



(4.6) 



Collecting the exponents of x in Eqs.(h5) and (f4 . 6|) we 
find for the two-qubit Hamiltonian: 



Ho = —J1S1S2 



h,Ai 



(4.7) 



This suggests that in general the magnetic field on spin 
i simulated by the algorithm takes the form: 



1 /V» A *+i' 

h< = A, I ° + 



2/3 



log 



Ji,Ai 



(4.8) 



whence the A^-spin Hamiltonian for an open ID chain 
becomes: 



N-l 



N-l 



H % = - ^ J i s i s i+i ~ Y h%Sl ~ ^NSN- (4.9) 



i=l 



i=l 



Note that hi can take any value for a given choice of fi- 
nite Ji and /3, by tuning the parameters Aj, A 4+ i. To 
prove that the algorithm sim ulates an Ising model with 
the Hamiltonian of Eq.(|4.9|) we proceed by induction. 
Assuming: 



i=N-l 
2 N -1 



L0 N 4>NX- H » [{s} * ] \{s}y) \ 



v=o 



L0 N 



1 w-l 

,0,Ai 1L L + 



-1/4 



(4.10) 



[where 4>n is defined in Eq.(2.14)] consider: 
1 



[Km lfl MHw> 



2 N -1 

2 N+1 -1 

L0N+1 Y <pN+lX- H ° N + ll{s} v ] \{s} y ), 
y=0 



The algorithm for the closed chain geometry takes a 
somewhat different form. Instead of applying i? Al one 
applies an ordinary ir/2 rotation on the first qubit, and 

closes the loop with , 1 . This results, as usual, in a 
superposition of ferro- and antifcrromagnetic last bonds, 
but also of positive and negative fields on si. To see this, 
consider: 



Sl -A 2S2 . IV) = II sf^R 1 \{-} N } = 



i=N-l 



N-l 



i=l 



where 



n {c^c^y^'y: ^x- H '^\{ S}y ) (4.12) 

v=o 



N-l 

h' n = - Ji 



SiS i+ i 



4=1 



N-l 



^ + (1/2/3) log(4- A Vc/- A2 ) - A 



NSn- 



When we now introduce a workbi t and apply f2, we find, 
after a calculation similar to Eq.(2.38): 



TV 

^' |Aii iV'>k=->=nK i ' As+i ' 

2 = 1 



JiAi+i 



-1/4 



2 N -1 

Y E <S>NX~ H »^\{ S }y)\w), 

w— ±1 y— 



(4.13) 



where A 



/v+i 



H 



N 



Ai, 

N 

i=i 



TV 



(4.14) 



i?^- is seen to be the correct Hamiltonian for a closed 
loop in the pres ence of the local fields hi. The sum over 
w = ±1 in Eq.(4.13) is such that w — sign Jjy =signAi, 
so that indeed, the algorithm produces a superposition 
over bond and field signs. Selecting a particula r sig n can 
be done with the interference method of Sec.IIF, and 



(4.11) 



where we used Eqs.(4.2) and ( |4.6| )-( |4.8| ). This proves the 
algorithm for the open chain case. 



the plaquette thus generated can be integrated into a 
higher dimensional lattice. A bond connecting plaquettes 
should not have to include a field term, since it presum- 
able connects spins which already have a field on them 
from the plaquette fabrication stage. Thus the situation 
in terms of controlling the introduction of a magnetic 
field is better than that of the bonds: arbitrary fields 
can be generated by the algorithm with full control over 
the field at every lattice point. It is interesting to point 
out in this context that it is known that the 2D fully 
antiferromagnetic Ising model with equal interactions, in 
the presence of a constant magnetic field is an NP-hard 
problem Eaj. 
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To conclude, we have introduced a new approach for 
the numerical study of statistical mechanics of Ising spin 
systems on quantum computers. The approach consists 
of an algorithm which allows one to construct a superpo- 
sition of qubit states, such that each state uniquely codes 
for a single configuration of Ising spins. Some stages of 
the algorithm, such as the construction of the open ID 
chain, are equivalent to the Markov process used in the 
transfer matrix formalism. Closing loops can be done by 
repeated measurements, requiring ~ exp(4/9J) attempts 
per frustrated plaquette. A crucial stage in the present 
algorithm, which leads to a polynomial increase in effi- 
ciency over algorithms which are based on measurements 
alone, is the use of an interference transformation. This 
transformation eliminates part of the superposition and 
thus determines whether the given plaquette will close in 
a frustrated or unfrustrated configuration. This is done 
in one step compared to ~ exp(4/3J) attempts per frus- 
trated plaquette in algorithms based on measurements 
alone. A central feature of the algorithm is that the 
quantum probability of each state in the superposition 
is exactly equal to the thermodynamic weight of the cor- 
responding configuration. When a measurement is per- 
formed, it causes the superposition to collapse into a 
single state. The probabilities of measuring states are 
ordered by the energies of the corresponding spin config- 
urations, with the ground state having the highest prob- 
ability. Therefore, statistical averages needed for calcu- 
lations of thermodynamic quantities obtained from the 
partition function, are approximated in the fastest con- 
verging order in the number of measurements. Unlike 
Monte Carlo simulations on a classical computer, consec- 
utive measurements on a quantum computer are totally 
uncorrelated. 

The algorithm applies to a large class of Ising systems, 
including partially frustrated models. A magnetic field 
can be incorporated as well without increase in the com- 
plexity, which is linear in the number of spins and bonds. 
The main problem of the algorithm is the limited con- 
trol it offers in the construction of a specific realization 
of bonds on the Ising lattice. An attempt to control all 
the bonds (and not only the prefabricated ones) by re- 
peating measurements may result in an exponential slow- 
down in performance as the temperature is lowered, and 
for this reason the algorithm fails to address the question 
of whether polynomial time (P) equals NP on a quantum 
computer, in the context of finding the spin glass ground 
state. In summary, this paper provides tools for the sim- 
ulation of Ising spin systems on a quantum computer, 
as efficiently as the best classical algorithms. Work em- 
ploying these tools to achieve speedup over classical al- 
gorithms is in progress. 
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FIGURES 
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FIG. 1. Scheme for numbering vertices on the binary Bethe tree, used in the Hamiltonian of Eq.(|3.1| 
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FIG. 2. All possible commutation relations of S (full arrows) and £1 (dashed arrows), (a) [Sij,Sik] = (needed for the 
Bethe lattice), (b) (I)-(VI) All combinations of Q, operators commute, (c) The commuting combinations of 5* and SI: (I) 
[Sij, Slikw] = 0, (II) [Sij,£lkiv>] = 0, (III) [Sij,£ljkw] = 0, (VI) [Sij,Qkiw] = 0. (d) The non-commuting combination of S and 
fl. (e) Additional commuting combinations needed in 4D and higher. 
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FIG. 3. a) Preparation of a single plaquette by three S operators (ferro (+), ferro, antiferro (-)) and closing with Q. This 
is followed by an interference transformation which erases the antiferromagnetic subspace, leaving in this case a frustrated (F) 
plaquette. b) A lattice with full density of prefabricated plaquettes (some unfrustrated - U), connected by f2 operators (dashed 
arrows), before measurement of the bonds. In this case all vertical bonds are present and the lattice is diluted in the horizontal 
bonds. 



